function y = func_ns(q, zeta, m, CONSTS)

    k0 = CONSTS.k0;
    a  = CONSTS.a;
    c  = CONSTS.c;

    y_vec = func_total(q, zeta, m, CONSTS);
    y_s_vec = - ((2*pi*k0^2*a)/c) * func_1(q, zeta, m, CONSTS) + ...
                       (2*pi/c*a) * func_2(q, zeta, m, CONSTS);
    y = y_vec - 1i * y_s_vec;

end